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ABSTRACT We describe an algorithm to 
compute native structures of proteins from 
their primary sequences. The novel aspects of 
this method are: 1) The hydrophobic potential 
was set to be proportional to the nonpolar 
solvent accessible surface. To make computa- 
tion feasible, we developed a new algorithm to 
compute the solvent accessible surface areas 
rapidly. 2) The supersecondary structures of 
each protein were predicted and used as re- 
straints during the conformation searching 
processes. This algorithm was applied to five 
proteins. The overall fold of these proteins can 
be computed from their sequences, with devia- 
tions from crystal structures of 1.48-4.48 A for 
C„ atoms. Proteins 31:247-257, 1998. 
©1998Wiley-Liss, Inc. 
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INTRODUCTION 

The success of protein structure prediction de- 
pends on our knowledge of protein structure, interac- 
tions, and folding mechanism. At present, there are 
only a few reliable conclusions: 1) Native structures 
of proteins are compact and have well-packed cores 
which are highly enriched in hydrophobic resi- 
dues. 1 - 2 2) Hydrophobic interaction is the driving 
force for protein folding. 3 - 4 Native structures of pro- 
teins have minimal solvent-exposed nonpolar sur- 
face areas. 5 3) Globular proteins are organized as a 
structural hierarchy, 6 - 7 i.e., secondary structure, su- 
persecondary structure, tertiary structure, and qua- 
ternary structure. 4) The proteins employ folding 
pathways to avoid extensively searching the whole 
conformation space. They fold by hierarchic conden- 
sation. 7 The folding pathway is suggested to be 
secondary structures, supersecondary structure, do- 
mains, and ultimately whole protein monomers. 8-10 

We developed a protein structure prediction algo- 
rithm based on this crude knowledge. First, the 
supersecondary structures were predicted with an 
artificial neural network method. 11 Then we searched 
for low-energy structures in the conformation space 



under the constraints suggested by the supersecond- 
ary structures. The energy function is very simple 
and has only two terms — a hydrophobic interaction 
term, and a van der Waals interaction term. We used 
a genetic algorithm to search the conformation space. 
According to our design of this model, hydrophobic 
interactions drive the peptide chain to fold; van der 
Waals forces are used to reject the incorrect compact 
structures during the hydrophobic collapse. Only the 
structures in which there is minimal conflict be- 
tween the hydrophobic interactions and the van der 
Waals interactions can survive and become domi- 
nant during the competition and selection process of 
the genetic algorithm. This algorithm was applied to 
five proteins. The overall fold of these proteins were 
computed from their sequences by this algorithm, 
with deviations from crystal structures of 1.48-4.48 
A for C a atoms (Table I). 

There have been several important advances in 
computer algorithms intended to predict native 3- 
dimensional structures of globular proteins from 
their amino acid sequences using simple energy 
functions. 12 ~ M The novel aspects of our method are: 
1) We set the hydrophobic potential to be propor- 
tional to the nonpolar solvent-accessible surface area 
(NSASA). Although this is a reasonable way of 
including the hydrophobic effects, it was not used in 
previous protein structure prediction algorithms. 12-14 
The computation of solvent-accessible surface area 
was very time-consuming. It would be prohibitively 
slow to incorporate it into protein structure predic- 
tion schemes which need to sample a large number of 
structures. We developed a new algorithm to calcu- 
late the solvent-accessible surface area rapidly. With 
this method we can compute the NSASA of every 
sampled structure in an acceptable time. 2) The 
supersecondary structures were predicted and used 
to derive soft constraints for the conformation search 
process. The identification of such structures repre- 
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TABLE I. Summary of the Computed Proteins 
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*The native disulphic 
simulation. 


e bond constraints were not used in the 



sents important progress along the folding pathway 
The conformation space is greatly reduced with these 
constraints. 

METHODS 
Review of Supersecondary 
Structure Prediction 

Supersecondary structure is defined as the combi- 
nation of two secondary structural elements with a 
short connecting peptide between one to five resi- 
dues in length. A short connecting peptide can have a 
large number of conformations. They play an impor- 
tant role in defining protein structures. A connecting 
peptide usually changes the trend of the protein 
backbones so as to form an antiparallel turn, a 
vertical corner, a twist, or just a slight bend in a 
peptide chain. 11 The conformations of the residues in 
the short connecting peptides are classified into five 
major types, namely, a, b, e, 1, or t, 16 each repre- 
sented by a region on the 4>-*(i map, respectively (see 
Fig. lb in Sun and Jiang 16 ). Supersecondary struc- 
tures are classified according to their component 
secondary structural elements, the length of the 
connecting peptide, and the type of residues in the 
connecting peptide. In a survey of 240 proteins, 16 it 
was found that there are 34 types of supersecondary 
structures which occur more than 5 times. Of these 
34 types there are 11 types of supersecondary struc- 
tures which occur more than 25 times. These 11 
types of supersecondary structures are called fre- 
quently occurring supersecondary structures. The 
34 types of supersecondary structures occurred alto- 
gether 766 times, among which the 11 frequently 
occurring supersecondary structures occurred 568 
times. This result shows that about 75% of the 
short-connecting peptides which occurred more than 
5 times belong to the 1 1 types of frequently occurring 
supersecondary structures. Sun et al. n developed an 
artificial neural network method to predict the 11 



frequently occurring supersecondary structure: 
H-b-H, H-t-H, H-bb-H, H-ll-E, E-aa-E, E-ea-E. 
H-lbb-H, H-lba-E. E-aal-E. E-aaal-E. and H-l-E. 
where "H H and "E" represent a-helix and p-strand, 
respectively. Each of these corresponds to a well- 
defined 3-dimensional motif (see Fig. 4 in Sun and 
Jiang 16 ). The method of Sun et al. was used for 
supersecondary structure prediction. The predicted 
supersecondary structure will not be rigidly imposed 
on the conformation. Rather, it will serve to define 
suitable constraints (Table II) on affected torsion 
angles. In this way, the size of the conformation 
space is greatly reduced. However, under these con- 
straints the conformation is still highly flexible and 
the structure can take on various shapes that are 
vastly different from the native shape. 

Peptide Chain Representation 

Amino acids are represented at the united-atom level. 
Bond lengths and bond angles are always fixed at their 
ideal values (according to Biosym's residue library). All 
the peptide bond dihedral angles are fixed in the trans 
(u> = 180°) conformation. The degrees of freedom in this 
reduced representation are the backbone and sidechain 
torsion angles <|>, 4*. and x (some residues have more 
than one sidechain torsion angle). 

Potential Energy Function 

Our potential function has two terms: a hydropho- 
bic interaction and a van der Waals interaction term, 

Etotal = &HH + Evdw* 

We define polar and nonpolar united atoms by 
their heavy atoms: carbon and sulphur are nonpolar; 
nitrogen and oxygen are polar. The hydrophobic 
potential is proportional to the solvent-accessible 
surface of nonpolar atoms, 

E HH = C h • NSASA 

where C h is a constant and is set to 0.031 and 
NSASA (in units of A 2 ) is the nonpolar solvent 
accessible surface area. 
We use a cut-off of 8 A for van der Waals interactions. 



where C v is a constant that is set to 0. 1 , ry is the 
distance between atom iand atom j, /?/and Rj are the 
van der Waals radii of atom i and atom j, and the 
summation is over all pairs of atoms with r f j < 8A. 
The function f vdtv is a van der Waals potential with a 
tapering-off at short distances (Fig. 3): 



2 

7* 



(r>0.8) 
(rs0.8). 
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Fig. 1 . The modified van der Waals potential energy function. 



This tapered van der Waals potential will not 
reject a structure with low hydrophobic energy for 
only a few steric conflicts. 

A backbone-to-backbone hydrogen bonding term is 
added for p-sheet structures. The distance between 
N -0 pair should be no more than 3.5 A and the 
out-of-plane dihedral angle between the oxygen and 
the peptide plane of the nitrogen (C-N-CJ should 
not exceed 40°. If these criteria are satisfied, an 
H-bond energy of -5.0 units is assigned. 

Two-Level Lattice Method for 

Solvent- Accessible Surface Area Calculation 

The hydrophobic effect is considered a principal 
force in the formation of the native protein struc- 
tures. 3 - 4 A reasonable way of including the hydropho- 
bic effect is to set the hydrophobic potential as 
proportional to the solvent-accessible surface area 
(SASA) of nonpolar atoms. Many algorithms for 
calculation of the SASA and molecule surface have 
been developed during the past 25 years. 17-40 Several 
fast numerical algorithms have been described in 
recent years. 38 39 However, calculation of the SASA 
has not been used in the folding simulation of whole 
protein molecules, in which a large number of SASA 
need to be calculated. In this work, we used a new 
algorithm to calculate the nonpolar solvent-acces- 
sible surface areas of about 300,000 conformations 
for each protein. Our method is based on the cube 
algorithm. 23-27 The whole molecule was put in a 
rectangular box. We introduced two levels of cubic 
lattice in the box, one was the coarse lattice (edge 
length 0.5 A), the other was the fine lattice (edge 
length 0. 1 A), by which each cube in the coarse lattice 
was divided into 125 sub-cubes. (The cube in the fine 
lattice is called a "sub-cube," which differs from the 
cube in the coarse lattice.) Each of the 125 sub-cubes 
was assigned a number from 0 to 124. Water mol- 
ecules were imitated by balls with a radius of 1.4 A. 
Before calculating the SASA, we built a library of the 



cubic decomposition of the accessible surfaces and 
inner parts of each kind of atoms. We selected the 
center of a particular cube to be the origin. Then we 
put a "hydrated" sphere with a radius equal to the 
sum of the van der Waals radii of the atom and a 
water molecule around the atom at the center of one 
of the sub-cubes in this cube. If the center of a cube 
was covered by the sphere, it was marked by "V" 
(V-cube). Then every V-cube was checked. If a V-cube 
was on the surface, this meant that at least one of its 
six neighboring cubes was not a V-cube, then it was 
marked by "S" (S-cube). If a V-cube was not on the 
surface, it was marked by "I" (I-cube). Thus, we have 
two kinds of V-cubes, i.e., S-cubes, which were on the 
surface of the "hydrated" sphere, and I-cubes, which 
were at the inner part of the "hydrated" sphere. The 
totality of S-cubes was an approximation of the 
surface of the "hydrated" sphere. The totality of 
I-cubes was the cubic decomposition of the inner part 
of the "hydrated" sphere. The positions (lattice coor- 
dinates relative to the origin) of the S-cubes and 
I-cubes were recorded in the library In such a way, 
the cubic decomposition of the "hydrated" sphere 
(include the surface and the inner part) whose center 
was at each of the 125 sub-cubes in the selected cube 
was recorded in the library according to the order of 
the sub-cube number. 

With this library, the SASA can be calculated 
rapidly: 

a) A protein molecule was put into the two-level 
lattice. For each atom, we determined which cube 
and sub-cube it was in. In other words, the 
Cartesian coordinates were transferred to lattice 
coordinates. 

b) If an atom was in the cube at (lx, ly, lz), which 
were the lattice coordinates of the center of the 
cube, and in the sub-cube whose number is n, 
then we look for the record of the cubic decomposi- 
tion of the "hydrated" sphere surface of this atom. 
The record was for the cube at the origin, so we 
translated them to (lx, ly, Iz). In such a way, we 
put the cubic approximation of the "hydrated" 
sphere surface of every atom in the lattice. The 
cubes that were occupied by the "hydrated" sphere 
surface were marked by "S." 

c) In the same way, we put the cubic approximation 
of the inner part of each atom on the lattice. These 
cubes were marked by "I." So, the S-cubes which 
were covered by the inner part of other "hydrated" 
spheres were remarked "L" 

d) Every S-cube was checked to ensure that it was 
really on the surface of the "hydrated" protein 
molecule. 

e) We counted the number of S-cubes. The number 
was proportional to SASA. Similarly, the total 
number of S-cubes belonging to nonpolar atoms 
would be proportional to NSASA. 
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(a) 1ROP 

1-50 1 2 3 4 5 

SEQUENCE MTKQEKTALNMARFIRSQTLTLLEKLNELDADEQADICESLHDHADELYR 

PREDICTED HHHHHHHHHHHHHHHHHHHHHHHHHHHH HHHHHHHHHHHHHHHHHHH 

NATIVE HHHHHHHHHHHHHHHHHHHHHHHHH HHHHHHHHHHHHHHHHHHHH 

51-56 

SEQUENCE SCLARF 

PREDICTED HHHHHH 

NATIVE HHHHHH 

(b) lUTG 

1-50 1 2 3 4 5 

SEQUENCE GICPRFAHVIENLLLGTPSSYETSLKEFEPDDTMKDAGMQMKKVLDSLPQ 

PREDICTED HHHHHHHHHHHHHHHbbHHHHHHHHHHHlbbHHHHHHHHHHHHHHHHbbH 

NATIVE HHHHHHHHHHH bbHHHHHHHHHH IbbHHHHHHHHHHHHHHHHbbH 

51-70 6 7 

SEQUENCE TTRENIMKLTEKIVKSPLCM 

PREDICTED HHHHHHHHHHHHHHHHHHH 

NATIVE HHHHHHHHHHHHHHH 

(C) 1CRN 

1-50 1 2 3 4 

SEQUENCE TTCCPSIVARSNFNVCRLPGTFEAICATYTGCIIIPGATCPGDYAN 

PREDICTED EEE HHHHHHHHHHHHHlbb HHHHHH HH1 EE EE HHHHH 

NATIVE EEEE HHHHHHHHHHHHH lbbHHHHHHHHl E E E E 

(d) 1R69 

1-50 1 2 3 4 5 

SEQUENCE SISSRVKSKRIQLGLNQAELAQKVGTTQQSIEQLENGKTKRPRFLPELAS 

PREDICTED HHHHHHHHHHHHHlbbHHHHHHHHlbbHHHHHHHHH HHHHHHHH 

NATIVE HHHHHHHHHHHH lbbHHHHHH lbbHHHHHHHHH HHHHHH 

51-63 6. . . 

SEQUENCE ALGVS VDWLLNGT 

PREDICTED HHlbbHHHHH 

NATIVE H lbbHHHHHH 

(e) 1CTF 

53-102 6 7 8 9 10. 

SEQUENCE E FDV I LKAAGANKVAV I KAVRGATGIiGLKE AKDLVE S APAALKEG VS KDD 

PREDICTED EEEEEEE HHHHHHHHHHHHlbbHHHHHHHHHHEEEEEEE HHH 

NATIVE EEEEEEE HHHHHHHHHHHHH lbbHHHHHHHHH EEEEE HHH 

103-120 11 12 

SEQUENCE AEALKKALEEAGAEVEVK 

PREDICTED HHHHHHHHHHH IEEE EEE 

NATIVE HHHHHHHHHHH 1 EEEEE 



Fig. 2. Comparison of the predicted supersecondary structures (predicted) and the X-ray- 
elucidated supersecondary structure (native) of (a) 1ROP, (b) 1UTG, (c) 1CRN, (d) 1R69, and (e) 
1CTF. 



The advantage of this method is that it does not 
require any calculation of distances between points. 
Besides table look-up. it only needs to perform about 
6N floating point multiplications in transferring the 
Cartesian coordinates to lattice coordinates, where 
N is the number of atoms. With this method we can 
calculate the SASA of N-terminal domain of the 434 
repressor (1R69, 484 nonhydrogen atoms) to within 



1 .2% accuracy in an average CPU time of 0. 1 1 5 s on a 
SGI PowerChallege R10000 processor at 194Mhz. 

Conformational Space Searching Algorithm 
and Supersecondary Structure Constraints 

We use genetic algorithm (GA) 42 - 43 to search the 
peptide chain conformational space for low-energy 
structures. In recent years there have been many 
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TABLE II. Corresponding Regions of the 
Supersecondary Structure Constraints? 



Supersecondary 
structures 


4> 




♦ 




a-helix 


-75° ~ 


-55° 


-50° ~ 


-30° 


p-strand 


-130° ~ 


-110° 


110° ~ 


130° 


a 


-150° ~ 


-30° 


-100° ~ 


50° 


b 


-230° ~ 


-30° 


100° ~ 


200° 


e 


30° ~ 


130° 


130°~ 


260° 


1 


30° ~ 


150° 


-60° ~ 


90° 


t 


-160° ~ 


-50° 


50° ~ 


100° 


undefined* 


-180° ~ 


0°** 


-180° ~ 


180° 



tThe rectangles were used as substitutes for the irregular 
regions of a. b. e, 1, and t on the <t>-i|> map. The most populated 
area of each region was Included in the corresponding rect- 
angle. 

*The residues that did not belong to any predicted supersecond- 
ary structures were classified as undefined. 
**For glycine this should be - 180° - 180*. 

studies of the use of GAs for protein structure 
prediction and other related structure optimization 
problems. 12 - 44-49 The basic idea of the genetic algo- 
rithm is to give better chances of survival and 
reproduction to the good individuals of the popula- 
tion (in our case, the low-energy structures). In this 
way the good genes (structural factors) will accumu- 
late and combine gradually to dominate the whole 
population. In the GA used in this work, a chromo- 
some consists of all the free variables in our peptide 
chain representation, i.e., it encodes the set of 4>, 
and x- Most residues have one, two, or three sidechain 
torsion angles, so the length of the chromosome is 
about 4N, where N is the number of residues. There 
are many versions of GA in its applications to protein 
folding simulation. Our G A procedure is: 

The initial population 

The initial population size was 500. These struc- 
tures were built by randomly selecting the backbone 
and sidechain torsion angles in the constrained 
regions. The backbone torsion angles 4> and of a 
residue were sampled uniformly in a certain region 
with a fineness of 1°. This region was defined by the 
position of the residue in the predicted super- 
secondary structures (Table II). The backbone tor- 
sion angles of the residues in the short connecting 
peptide of any super-secondary structure were con- 
strained to lie in the corresponding regions on the 
<t>-*|> map (Fig. 1). The backbone dihedral angles of 
the residues in the a-helix and p-strand were re- 
strained to within ±10* of their ideal value, which 
was set to (-65°, -40°) and (-120°, 120°) respec- 
tively. In this way. predicted supersecondary struc- 
tures were used to greatly reduce the allowable 
variation of the affected torsion angles. For the other 
residues not predicted to lie in any supersecondary 
structure, their backbone torsion angles are sampled 
randomly on the left half (<(> < 0) of the 4>-i|i map (for 



glycine, it was the whole <M* map). For the sidechain 
torsion angles, the constraints were based on the 
sidechain rotamer library. 41 The mean value of the 
sidechain torsion angle in the rotamer library was 
selected according to its occurring ratio. After a 
mean value was selected, an integer value is selected 
randomly from the interval [mean value — standard 
deviation, mean value + standard deviation]. These 
500 structures were the parent individuals of the 
first generation. 

Fitness criterion 

The potential energy of the 500 parent individuals 
were computed. Potential energy was used as the 
objective function. Then we mapped the objective 
function onto a fitness scale. If there are a few 
extraordinary individuals in the early stage of the 
GA process, they will take over a significant propor- 
tion of the population after several generations. This 
is a leading cause of premature convergence. On the 
other hand, if there is still significant diversity 
within the population in the later stages of the GA 
process, then the population average fitness may be 
close to the population best fitness and the best 
members may not become dominant in the popula- 
tion. In this case, the GA process becomes a random 
walk. 43 To prevent premature convergence and ran- 
dom walk, we used a generation-dependent fitness 
scaling. 

Fitness gnJ = 1 + C gn -= — = 

gn. max *-'gn,m\n 

C g n = C 0 + incr * gn 

where Fitness gn j is the fitness of the ith individual in 
gnth generation, E glhmax , E g „ tmJnt and E gnJ is highest, 
lowest, and the ith individual's potential energy in 
the gnth generation, Cq is a constant that is set to be 
0.02, incr is increment of the ratio of fitness of the 
best individual (with lowest energy) to the worst 
individual (with highest energy in each generation). 
We set incr= 0.0016 in our computations. 

In each generation, the fitness of the worst indi- 
vidual was always set to 1 , the best individual was 
1 + Cb + incr • gn. This scaling strategy is a variant 
of the fitness scaling methods that focus on the ratio 
of the fitness of the best individual and the average 
fitness. Premature convergence and random walk 
can be prevented by this scaling strategy. 

The crossover operation 

Pairs of individuals were selected randomly for 
crossover operation. The probability for an indi- 
vidual to be selected was f/lf, where f f was the 
fitness of the individual and 2/Vas the summation of 
the fitness of all the individuals in the population. 
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NATIVE PREDICTED 

Fig. 3. Structure comparison between the crystal structure (native) and the computed structure 
(predicted) of repressor of primer. 



The chromosome looked like: 

<t> 1 + 1 X 11 X 1 2l4>2^2X2 1 X22X23I 4>3^3X3 1 1 

I 4>ntnXnlXn2l I+N+NXNiXN* 

where n was the residue number, Xni and Xn2 was the 
Xi and X2 of the nth residues, | is a possible crossover 
site. In order to keep the correlation between <t> and *|/ 
in the super-secondary structures, and the correla- 
tion of xi. X2. and X3 in the sidechain rotamer library, 
the crossover site is not allowed to occur between the 
4>, *Ji, and x of the same residue. A selected pair of 
chromosomes would undergo a fixed number (chosen 
to be 1 in this study) of crossovers at randomly 
chosen allowable sites. 

The mutation operation 

Two kinds of mutation operators were used. The 
first mutation operator may change the conforma- 
tion dramatically. When this operator acted on a 
peptide chain, all the values of the backbone and 
sidechain torsion angles of a randomly chosen resi- 
due were reselected from their corresponding con- 
strained regions. We made a copy of the 500 parent 
individuals and modified this copied population Mj 
times, each time by applying the operator to a 
randomly selected individual from this population. 
An individual can be selected more than one time, so 
there may be changes in torsion angles in more than 
one residue in a chromosome. The second mutation 
operator is for a more local search of conformational 
space. 12 It will perturb some residues' torsion angles 



(4>, and x) by a random angle between —5° and 5°. The 
number of perturbed residues of each individual is A/ 2 . 
This operator was also applied to every parent indi- 
vidual so that in total 500 offspring were produced. 

Again, for the purposes of preventing premature 
convergence and random walk, we made M x and M 2 
decrease as the search proceeds: 

M, = 1 + P- exp(-gn/Neff) 

AA 

M z = 1 + — ■ exp(-gn/Neff) 

where Pis set to 500, Af is the number of residues in 
the protein, gn is the generation, and Neff is a 
constant set to 150. 

Selection 

Now the population consists of 500 parent conforma- 
tions, 500 crossed offspring, and 1,000 mutated off- 
spring. The total population is 2,000. The potential 
energy of these 2,000 conformations were computed and 
only the 500 lowest-energy conformations were selected 
into the next generation as parent conformations. 

Convergence 

At least 100 generations of GA were performed for 
each protein. After 100 generations, the GA process 
will stop only if the decrease of the lowest energy in 
the population is less than 1 unit during the last 20 
generations. On average, about 150 generations of 
GA were performed for each protein. 
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RESULTS 
Repressor of Primer (1ROP) 

Repressor of primer is a 4-helix bundle protein 
that is composed of two identical monomers. Each 
monomer has 56 residues and forms a a-turn-a 
structure, which does not belong to the 1 1 frequently 
occurring supersecondary structures. The predicted 
secondary structures (Fig, 2a) were used as con- 
straints. After computing the conformation using our 
algorithm, we calculated its distance matrix error 
(DME) to the crystal structure. The computed struc- 
ture matches the crystal structure with a DME of 
1.48 A (Fig. 3). 

Uteroglobin (1UTG) 

Uteroglobin is a 4-helix protein that has 70 resi- 
dues. The predicted supersecondary structures are 
tx-bb-a-lbb-a-bb-a (Fig. 2b), With these supersecond- 
ary structure constraints we computed the 3- 
dimensional structure. The computed structure 
matches the crystal structure with a DME of 3.47 A 
(Fig. 4). 

Crambin (1CRN) 

Crambin is a 46-residue protein with two a-helix 
and a pair of p-strands. It has three disulphide 
bonds. We did not use the disulphide bond con- 
straints. The predicted supersecondary structures 
are p-loop-a-lbb-ct-l-(J-loop-a (Fig. 2c). The computed 
structure matches the crystal structure with a DME 
of2.73A(Fig. 5). 



N-Terminal Domain of the 434 
Repressor (1R69) 

The crystal structure of the N-terminal domain of 
the 434 repressor has 63 residues and is composed of 
five helices. The predicted supersecondary struc- 
tures are tx-lbb-ce-lbb-a-loop-a-lbb-a (Fig. 2d). The 
computed structure matches the crystal structure 
with a DME of 4.48 A (Fig. 6). 

C-Terminal Domain of the L7 (SLASH) *L 12 50 
S Ribosomal Protein (1CTF) 

This protein has 68 residues. It has six secondary 
structures — three a-helix and three p-strands. This 
protein is the most complex example in this study. 
The predicted supersecondary structures are p-a-lbb- 
a-p-a-l-p (Fig. 2e). The computed structure matches 
the crystal structure with a DME of 4.00 A (Fig. 7). 

DISCUSSION 

The predicted super-secondary structures and the 
native supersecondary structures of these five pro- 
teins are shown in Fig. 2. In these five proteins, there 
are 21 secondary structures and 16 short connecting 
peptides. Ten short connecting peptides were identi- 
fied to be in one of the 11 frequently occurring 
supersecondary structures. Most of the supersecond- 
ary structures are correctly predicted. For these five 
proteins the correctness ratio is 90.1%. Although the 
accuracy is high, in some instances the predicted 
structures do not align precisely with those observed 
in the crystal structures. If the backbone torsion 
angle (<J>, i|>) of a few consecutive residues were 
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NATIVE PREDICTED 

Fig. 5. Structure comparison between the crystal structure (native) and the computed structure 
(predicted) of repressor of Crambin. 



NATIVE PREDICTED 

Fig. 6. Structure comparison between the crystal structure (native) and the computed structure 
(predicted) of N-terminal domain of the 434 repressor. 



restrained in wrong regions, the peptide chain may 
have a wrong trend at this segment. If the segment 
was at the central part of the peptide chain, the 
overall fold may be misdetermined. In this study, 
such a fatal mistake has not occurred in the supersec- 
ondary structure prediction. The most serious struc- 
ture distortion caused by the errors of the supersec- 
ondary structure prediction was in the crambin. At 
the C-terminal of the peptide chain, an incorrectly 
predicted a-helix (from residue 41 to 45) was im- 
posed on the peptide chain as a constraint (Fig. 4). 



As a result, a wrong structure was formed in this 
terminal (Fig. 7). 

In Figure 2 one can find that in some cases 
supersecondary structure was not correctly pre- 
dicted at only one or two residues, while the neighbor- 
ing residues were all restrained in the correct re- 
gions. In these cases, the residue the peptide chain 
will turn to a wrong direction at this point. But if the 
native-like structures are favored by the potential, 
the nearby residues will move to compensate for this 
mistake. As a result, a native-like profile can still be 
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NATIVE PREDICTED 

Fig. 7. Structure comparison between the crystal structure (native) and the computed structure 
(predicted) of C-terminal domain of the ribosomal protein L7/*L12. 



formed. For example, the computed structure of lutg 
(Fig. 4) was native-like, while at two central residues 
(15 and 28) the supersecondary structures were 
predicted incorrectly. This is possible because of the 
flexibility of the model. In this model, the peptide 
chains are more flexible than those in the fixed 
secondary structure model 12 where the affected tor- 
sion angles are fixed at their ideal values (without 
any degree of freedom). In contrast, the torsion 
angles in our model can rotate in some regions, 
which is determined by the predicted supersecond- 
ary structures and the sidechain rotamer library. 

With the supersecondary structure constraints the 
conformation space of the peptide chain is greatly 
reduced, but the DME of the peptide chain can still 
be large (over 15 A). For example, in the simulation 
of the 434 repressor we observed many conforma- 
tions with a DME exceeding 15 A, especially in the 
early generations. This indicates that the overall 
fold of a protein molecule cannot be well defined 
only by the supersecondary structures. Genetic algo- 
rithm was used to search the reduced conformation 
space for low-energy structures. Four of the five 
computed structures are similar to the correspond- 
ing X-ray elucidated structures. The DME of repres- 
sor of primer, uteroglobin, crambin, and C-terminal 
domain of L7(SLASH)*L12 50 S ribosomal protein 
are all smaller than or equal to 4.0 A. For the 
N-terminal domain of the 434 repressor, there are 
three supersecondary structures. One of the connect- 
ing peptides, a six-residue loop, cannot be recognized 
as belonging to any of the 1 1 kinds of supersecondary 
structures. The peptide chain is divided into two 
fragments which are connected by the loop region. 
The computed structure of each of these two frag- 



ments is similar to their corresponding parts in the 
crystal structure, but the relative position of the two 
fragments is incorrectly determined by the loop 
region. 

An ideal potential should give higher values to 
non-native conformations than the native conforma- 
tion. Recent studies 50-52 indicated that some poten- 
tials can distinguish between correct and certain 
incorrect structures with a high degree of success. 
However, in the protein folding simulation there is 
an astronomical number of candidates in the confor- 
mation space. If a potential can identify 99.99% 
non-native structures, while it gives 0.01% of them 
lower energy than the native structure, then there 
are still uncountably many minima with lower en- 
ergy than that of the native structure on the energy 
landscape. In this situation, our hope of finding a 
native-like structure lies in the possibility that most 
of the low-energy structures are "near" the native 
structure to form a cluster of native-like structures. 
This appears to be the case for the repressor of 
primer and uteroglobin, where the energy of the 
computed structures is lower than that of their 
native structures. The other three proteins are more 
complex; the energy of the computed structures are 
higher than that of their native structures. This is 
mainly caused by the van der Waals term and it 
indicates that a more efficient local search method is 
needed. 

In recent years, important progress has been made 
in computing the 3-dimensional structures of pro- 
teins from their sequences using simple energy 
function. The attractive aspect of this method is that 
a simple model should be much easier to improve 
than highly parameterized ones, and the prediction 
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results of these simple models is arguably compa- 
rable to those of the more complex models. 14 Encour- 
aging results have been reported by Sun et al. 12 They 
developed a model that predicted reasonably well 
the known tertiary folds of 7 out of 1 0 small proteins. 
Their method used experimental secondary struc- 
tures, in which the backbone dihedral angles (o>, 
are fixed at the ideal values. Three of these ten 
proteins are also considered in this study. They are 
repressor of primer, crambin, and N- terminal do- 
main of the 434 repressor. The DME reported in Sun 
et al. are 1.65 A, 4.87 A, and 5.55 A, respectively. In 
our results, the DME of these proteins are 1.48 A, 
2.73 A, and 4.48 A, respectively. This suggests that 
supersecondary structure constraints and better mod- 
eling of the hydrophobic interaction are of consider- 
able utility in protein structure computation. 

CONCLUSION 

One important step toward building a tertiary 
structure is to identify how secondary structures as 
building blocks arrange themselves in space. Good 
supersecondary structure prediction methods can 
provide important information in the prediction of 
protein tertiary structure. 

The structure of a protein is determined by the 
competition and cooperation of all of the interac- 
tions, especially hydrophobic interaction and van der 
Waals interaction. Correctly including the hydropho- 
bic interaction is extremely important. Although the 
nature of hydrophobic interaction is not -completely 
understood, it is suggested that protein-solvent inter- 
action depends on the solvent-accessible surface 
area of the protein molecule. 53,54 

The goal of this study is to suggest a way to 
capture these two main features in our current 
understanding of protein structure, interaction, and 
folding mechanism. The results show that some 
small protein structures can be determined by a 
model that carefully adduces these points. 
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